HW
1. Problem 11.1
data <- read_excel("Chapter11_exercises_data.xls")## New names:
## • `` -> `...8`
## • `` -> `...9`data$GSF <- as.numeric(data$GSF)
data$GSJ <- as.numeric(data$GSJ)
data <- data %>%
dplyr::select(Date, GSF, GSJ) %>%
na.omit()# Split into estimation and prediction samples
train <- data[1:(nrow(data) - 20), ] %>% dplyr::select(GSF, GSJ)
test <- data[(nrow(data) - 19):nrow(data), ] %>% dplyr::select(GSF, GSJ)
summary(train)## GSF GSJ
## Min. :-4.6800 Min. :-3.232
## 1st Qu.: 0.2828 1st Qu.: 0.400
## Median : 2.2296 Median : 1.635
## Mean : 2.1558 Mean : 1.892
## 3rd Qu.: 3.8616 3rd Qu.: 3.311
## Max. :10.5140 Max. : 9.799summary(test)## GSF GSJ
## Min. :-8.7830 Min. :-8.760
## 1st Qu.:-4.4226 1st Qu.:-3.623
## Median :-0.8497 Median :-1.162
## Mean :-1.7085 Mean :-1.133
## 3rd Qu.: 0.5783 3rd Qu.: 1.406
## Max. : 6.4845 Max. : 4.582GSF and GSJ represent the quarterly house price growth rate calculated based on the original house price indices for San Francisco-Oakland-Fremont and San Jose-Sunnyvale-Santa Clara, respectively. For the train sample, GSF has a median of 2.2296 and mean of 2.1558. The middle 50% of data lies between 0.2828 and 3.8616. For the train sample of GSJ, the median is 1.635 and mean is 1.892, while the middle 50% of data lies between 0.400 and 3.311. In both test samples, the median and mean are both negative, and the min, max, and quartiles are all much lower compared to the train data. This is because of the dynamics of the time period from 2008 to 2012.
# Create TS objects using estimation sample
GSF_train <- ts(train$GSF,
start = c(1975, 2),
end = c(2007, 3),
frequency = 4)
GSJ_train <- ts(train$GSJ,
start = c(1975, 2),
end = c(2007, 3),
frequency = 4)
# VAR model
VAR_data <- window(ts.union(GSF_train, GSJ_train), start = c(1975, 2), end = c(2007, 3))
VARselect(VAR_data, lag.max = 10)## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 3 3 3 3
##
## $criteria
## 1 2 3 4 5 6 7 8
## AIC(n) 1.260763 1.259156 1.003914 1.028296 1.034471 1.071311 1.014798 1.077252
## HQ(n) 1.317364 1.353491 1.135982 1.198098 1.242007 1.316580 1.297802 1.397989
## SC(n) 1.400138 1.491447 1.329121 1.446420 1.545511 1.675267 1.711671 1.867041
## FPE(n) 3.528188 3.522788 2.729665 2.797876 2.816523 2.924191 2.766069 2.947888
## 9 10
## AIC(n) 1.071033 1.080588
## HQ(n) 1.429504 1.476793
## SC(n) 1.953739 2.056210
## FPE(n) 2.934116 2.967940We will use a lag length of 3 for the VAR model.
var_model <- VAR(y = VAR_data, p = 3)
summary(var_model)##
## VAR Estimation Results:
## =========================
## Endogenous variables: GSF_train, GSJ_train
## Deterministic variables: const
## Sample size: 127
## Log Likelihood: -418.202
## Roots of the characteristic polynomial:
## 0.9124 0.6813 0.6813 0.5674 0.319 0.319
## Call:
## VAR(y = VAR_data, p = 3)
##
##
## Estimation results for equation GSF_train:
## ==========================================
## GSF_train = GSF_train.l1 + GSJ_train.l1 + GSF_train.l2 + GSJ_train.l2 + GSF_train.l3 + GSJ_train.l3 + const
##
## Estimate Std. Error t value Pr(>|t|)
## GSF_train.l1 0.27804 0.12489 2.226 0.0279 *
## GSJ_train.l1 0.65073 0.14912 4.364 2.72e-05 ***
## GSF_train.l2 -0.16813 0.13269 -1.267 0.2076
## GSJ_train.l2 -0.33706 0.17370 -1.940 0.0547 .
## GSF_train.l3 0.57698 0.12499 4.616 9.88e-06 ***
## GSJ_train.l3 -0.08235 0.15007 -0.549 0.5842
## const 0.16928 0.21193 0.799 0.4260
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 1.653 on 120 degrees of freedom
## Multiple R-Squared: 0.6244, Adjusted R-squared: 0.6056
## F-statistic: 33.25 on 6 and 120 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation GSJ_train:
## ==========================================
## GSJ_train = GSF_train.l1 + GSJ_train.l1 + GSF_train.l2 + GSJ_train.l2 + GSF_train.l3 + GSJ_train.l3 + const
##
## Estimate Std. Error t value Pr(>|t|)
## GSF_train.l1 -0.06316 0.10844 -0.583 0.561319
## GSJ_train.l1 1.00212 0.12947 7.740 3.48e-12 ***
## GSF_train.l2 -0.19186 0.11521 -1.665 0.098453 .
## GSJ_train.l2 -0.29834 0.15081 -1.978 0.050200 .
## GSF_train.l3 0.37086 0.10852 3.417 0.000864 ***
## GSJ_train.l3 -0.01054 0.13030 -0.081 0.935652
## const 0.28633 0.18401 1.556 0.122327
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 1.435 on 120 degrees of freedom
## Multiple R-Squared: 0.637, Adjusted R-squared: 0.6189
## F-statistic: 35.1 on 6 and 120 DF, p-value: < 2.2e-16
##
##
##
## Covariance matrix of residuals:
## GSF_train GSJ_train
## GSF_train 2.733 1.688
## GSJ_train 1.688 2.060
##
## Correlation matrix of residuals:
## GSF_train GSJ_train
## GSF_train 1.0000 0.7112
## GSJ_train 0.7112 1.0000In the equation for GSF, the second lag of GSF and the third lag of GSJ have p-values higher than 0.05, so we could consider removing them from the model. The equation for GSF has an adjusted R-squared value of 0.6056, which is pretty good. In the equation for GSJ, the first and second lags for GSF and the third lag of GSJ have p-values higher than 0.05, so they may not have a meaningful influence on GSJ. The model has an adjusted R-squared value of 0.6189, similar to the first equation.
2. Problem 11.2
grangertest(GSF_train, GSJ_train, order = 3)## Granger causality test
##
## Model 1: GSJ_train ~ Lags(GSJ_train, 1:3) + Lags(GSF_train, 1:3)
## Model 2: GSJ_train ~ Lags(GSJ_train, 1:3)
## Res.Df Df F Pr(>F)
## 1 120
## 2 123 -3 4.0141 0.00921 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1grangertest(GSJ_train, GSF_train, order = 3)## Granger causality test
##
## Model 1: GSF_train ~ Lags(GSF_train, 1:3) + Lags(GSJ_train, 1:3)
## Model 2: GSF_train ~ Lags(GSF_train, 1:3)
## Res.Df Df F Pr(>F)
## 1 120
## 2 123 -3 6.5504 0.0003863 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Using an order of 3, the granger-causality tests both have a p-value below 0.05, so they suggest that GSF granger-causes GSJ and that GSJ also granger-causes GSF. This indicates that past values of GSF can be used to predict future values of GSJ, and vice versa.
3. Problem 11.3
plot(irf(var_model))
The first plot shows that a shock in GSF causes an immediate significant response in GSF followed by a gradual decline around the second lag. This suggests that the economy responds significantly to its own shocks in the short term and maintains a positive response throughout. The second plot similarly shows that a shock in GSF causes a positive and significant response in GSJ at first followed by a gradual decline to near baseline level around the second lag, indicating a connection between the two economies. The third plot shows that a shock in GSJ causes a positive response in GSF in the short term followed by a gradual decline to baseline by the third lag, at which point the response is about zero. The final plot shows that a shock in GSJ causes a positive and significant response in GSJ in the short term followed by a gradual decline to baseline by the third lag, at which point the response is roughly zero throughout.
4) Problem 7.8
a)
retaildata <- read_excel("retail.xlsx", skip=1)
myts <- ts(retaildata[,"A3349873A"],
frequency=12, start=c(1982,4))
autoplot(myts)
Multiplicative seasonality is necessary because the seasonal
fluctuations vary with time. The amplitude of this data grows with time,
signaling multiplicative tendencies.
b)
fit1 <- hw(myts, seasonal = "multiplicative")
fit2 <- hw(myts, seasonal = "multiplicative", damped = TRUE)
autoplot(myts) +
autolayer(fit1, series="HW multiplicative forecasts", PI=FALSE) +
autolayer(fit2, series="HW multiplicative forecasts - damped", PI=FALSE)
c)
rmse1 <- accuracy(fit1)["Training set", "RMSE"]
rmse1## [1] 13.29378rmse2 <- accuracy(fit2)["Training set", "RMSE"]
rmse2## [1] 13.30494The RMSEs are basically the same, but the one for the multiplicative method is slightly lower so I will prefer that one over the damped method.
d)
tsdisplay(fit1$res)
lbtest <- Box.test(fit1$res, lag=20, type="Ljung-Box")
lbtest##
## Box-Ljung test
##
## data: fit1$res
## X-squared = 35.369, df = 20, p-value = 0.01822No, the residuals from the best method do not exhibit white noise. There are a few lags that cross the bands, signalling autocorrelation. Also, the Ljung-Box test (p = 0.01822) confirms that the residuals are not white noise.
e)
myts.train <- window(myts, end=c(2010,12))
myts.test <- window(myts, start=2011)
fit3 <- hw(myts.train, h=36, seasonal = "multiplicative")
fc <- snaive(myts.train, h=36)
autoplot(myts) +
autolayer(fit3, series="multiplicative", PI=FALSE) +
autolayer(fc, series="snaive", PI=FALSE)
accuracy(fit3, myts.test)## ME RMSE MAE MPE MAPE MASE
## Training set 0.03021223 9.107356 6.553533 0.001995484 3.293399 0.4107058
## Test set 78.34068365 94.806617 78.340684 19.945024968 19.945025 4.9095618
## ACF1 Theil's U
## Training set 0.02752875 NA
## Test set 0.52802701 1.613903accuracy(fc, myts.test)## ME RMSE MAE MPE MAPE MASE
## Training set 7.772973 20.24576 15.95676 4.702754 8.109777 1.000000
## Test set 81.744444 100.00869 82.06667 20.549055 20.669738 5.143067
## ACF1 Theil's U
## Training set 0.7385090 NA
## Test set 0.6830879 1.67023rmse3 <- accuracy(fit3, myts.test)["Test set", "RMSE"]
rmse3## [1] 94.80662rmse4 <- accuracy(fc, myts.test)["Test set", "RMSE"]
rmse4## [1] 100.0087Yes, the holt-winters multiplicative method beats the seasonal naive approach. The rmse is lower and the graph shows it is slightly closer to the actual test set values than the snaive approach.
5) Problem 7.9
#Box-Cox transformation
lambda <- BoxCox.lambda(myts)
myts.bc <- BoxCox(myts, lambda)
myts.bc <- ts(myts.bc[,1], frequency=12, start=c(1982,4))
#STL decomposition
stl_decomp <- stl(myts.bc, s.window="periodic", robust = TRUE)
autoplot(stl_decomp)
#get seasonally adjusted data
seasonal_component <- stl_decomp$time.series[, "seasonal"]
seasonally_adjusted <- myts.bc - seasonal_component
#define training/test sets
myts_train <- window(myts.bc, end=c(2010,12))
myts_test <- window(myts.bc, start=2011)
#perform ETS
ets.fit <- ets(myts_train)
ets.forecast <- forecast(ets.fit, h = 36)
#check RMSE
rmse5 <- accuracy(ets.forecast, myts_test)["Test set", "RMSE"]
rmse5## [1] 0.5734974The RMSE of 0.5734974 is the lowest compared to the other best previous forecasts.
Question 6 11-a)
data(visitors)autoplot(visitors)
The data exhibits a multiplicative increasing trend. Additionally, a
seasonal pattern is evident, recurring every year.
11-b)
fc <- hw(subset(visitors,end=length(visitors)-24),
damped = TRUE, seasonal="multiplicative", h=24)
autoplot(visitors) +
autolayer(fc, series="HW multi damped", PI=FALSE)+
guides(colour=guide_legend(title="Yearly forecasts"))
11-c)
Due to the increase in the magnitude of the visitor numbers, a multiplicative seasonal model is necessary.
11-d)
visitors_train <- subset(visitors, end = length(visitors) - 24)
visitors_test <- subset(visitors, start = length(visitors) - 24)
# 1. ETS model
fc_ets <- forecast(ets(visitors_train), h = 24)
plot(fc_ets)
# 2. Additive ETS on Box-Cox transformed series
lambda <- BoxCox.lambda(visitors_train)
fc_ets_add_BoxCox <- forecast(ets(visitors_train, lambda = lambda, additive.only = TRUE), h = 24)
plot(fc_ets_add_BoxCox)
# 3. Seasonal naive method
fc_snaive <- snaive(visitors_train, h = 24)
plot(fc_snaive)
# 4. STL decomposition on Box-Cox transformed data with ETS on seasonally adjusted data
visitors_train_bc <- BoxCox(visitors_train, lambda)
stl_out <- stl(visitors_train_bc, s.window = 12) # Use seasonal argument
visitors_sa <- stl_out$time.series
plot(visitors_sa)
11-e)
checkresiduals(fc_ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(M,Ad,M)
## Q* = 20.677, df = 24, p-value = 0.6577
##
## Model df: 0. Total lags used: 24checkresiduals(fc_ets_add_BoxCox)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,Ad,A)
## Q* = 15.206, df = 24, p-value = 0.9146
##
## Model df: 0. Total lags used: 24checkresiduals(fc_snaive)
##
## Ljung-Box test
##
## data: Residuals from Seasonal naive method
## Q* = 295.02, df = 24, p-value < 2.2e-16
##
## Model df: 0. Total lags used: 24The Ljung-Box test reveals that the seasonal naive method fails to capture all the trends in the data. This is because its residuals exhibit a strong autocorrelation, indicated by a very high Q statistic and a very low p-value from the test. On the other hand, the Additive ETS model on Box-Cox transformed series emerges as the best candidate. This model demonstrates the lowest Q statistic and the highest p-value in the Ljung-Box test, suggesting minimal autocorrelation in the residuals and a good fit for the data.
11-f)
# first, make functions to make model to yield forecast class object
fets_add_BoxCox <- function(y, h) {
forecast(ets(
y,
lambda = BoxCox.lambda(y),
additive.only = TRUE
),
h = h)
}
fstlm <- function(y, h) {
forecast(stlm(
y,
lambda = BoxCox.lambda(y),
s.window = frequency(y) + 1,
robust = TRUE,
method = "ets"
),
h = h)
}
fets <- function(y, h) {
forecast(ets(y),
h = h)
}
# I'll compare the models using RMSE
sqrt(mean(tsCV(visitors, snaive, h = 1)^2, na.rm = TRUE))## [1] 32.56941sqrt(mean(tsCV(visitors, fets_add_BoxCox, h = 1)^2,
na.rm = TRUE))## [1] 18.8505sqrt(mean(tsCV(visitors, fstlm, h = 1)^2,
na.rm = TRUE))## [1] 17.49642sqrt(mean(tsCV(visitors, fets, h = 1)^2, na.rm = TRUE))## [1] 18.52985sqrt(mean(tsCV(visitors, hw, h = 1,
seasonal = "multiplicative")^2,
na.rm = TRUE))## [1] 19.62107tsCV errors show that the best model is the STL + ETS(M, A, N) model and the worst model is seasonal naive model. If I hadn’t calculated accuracy using test set, I couldn’t have known that the forecasts from seasonal naive method were the most accurate ones.
Question 7 13-a)
data(hsales)
is.ts(hsales)## [1] TRUEThe data transformation is necessary due to its non-stationary nature (lack of mean reversion).
13-b)
tsdisplay(hsales)
hsales_diff <- diff(hsales, lag = 12 , differences = 1)
tsdisplay(hsales_diff)
Using seasonal difference would take the seasonal pattern present in the data that repeat every 12 months into account.
13-c)
#hsales_diff
hsales_model <- arima(hsales_diff, order = c(13, 0, 0))
AIC(hsales_model)## [1] 1641.46hsales_model_2 <- arima(hsales_diff, order = c(25, 0, 0))
AIC(hsales_model_2)## [1] 1641.232Based on the exponentially decaying ACF with significant spikes at lags 13 and 25, I’ve chosen to evaluate ARIMA(13,0,0) and ARIMA(25,0,0) models. Since a lower AIC value indicates a better fit, ARIMA(25,0,0) would be a better model.
13-d
coef(hsales_model_2)## ar1 ar2 ar3 ar4 ar5 ar6
## 0.836300050 -0.018559397 0.007202887 0.004160205 0.068567763 0.104486076
## ar7 ar8 ar9 ar10 ar11 ar12
## -0.085498706 0.041765085 -0.079485039 0.093966567 -0.002479969 -0.581887389
## ar13 ar14 ar15 ar16 ar17 ar18
## 0.491495817 0.007722343 -0.014495612 0.024323854 0.061152752 0.006687802
## ar19 ar20 ar21 ar22 ar23 ar24
## -0.072474289 -0.054275164 0.008676560 -0.011662305 0.011350102 -0.259985765
## ar25 intercept
## 0.273731526 -0.189719054checkresiduals(hsales_model_2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(25,0,0) with non-zero mean
## Q* = 14.475, df = 3, p-value = 0.002325
##
## Model df: 25. Total lags used: 28The residuals exhibit an approximately normal distribution centered around zero, suggesting characteristics close to a white noise process.
13-e)
plot(forecast(auto.arima(hsales),h=24))
13-f)
hsales %>% ets() %>% forecast(h=24) %>% autoplot()
The ETS forecast exhibits a stronger seasonal pattern compared to the other models.
Question 8: 8.18.
a.
gdp.NA <- Quandl("FRED/GDP", api_key= "mVHzLbBisoncyRnYC27C", type = "ts", transform = "rdiff")
str(gdp.NA)## Time-Series [1:299] from 1947 to 2022: 0.0115 0.0147 0.0407 0.0231 0.0257 ...head(gdp.NA)## Qtr1 Qtr2 Qtr3 Qtr4
## 1947 0.01153131 0.01470516 0.04070757
## 1948 0.02308803 0.02568281 0.02432063b. Plot graphs of the data, and try to identify an appropriate ARIMA model.
autoplot(gdp.NA)
ndiffs(gdp.NA)## [1] 1ggtsdisplay(diff(gdp.NA))
auto.arima(gdp.NA)## Series: gdp.NA
## ARIMA(1,1,3)(2,0,0)[4]
##
## Coefficients:
## ar1 ma1 ma2 ma3 sar1 sar2
## 0.1792 -0.9982 0.2036 -0.1514 -0.0820 -0.1125
## s.e. 0.2982 0.2963 0.2601 0.0632 0.0639 0.0814
##
## sigma^2 = 0.0001602: log likelihood = 881.22
## AIC=-1748.44 AICc=-1748.05 BIC=-1722.56From auto.arima, ARIMA(1,1,3)(2,0,0) model will fit well to the data.
c. Do residual diagnostic checking of your ARIMA model. Are the residuals white noise?
gdp.NA_arima <- auto.arima(gdp.NA)
autoplot(gdp.NA_arima$residuals)
checkresiduals(gdp.NA_arima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,3)(2,0,0)[4]
## Q* = 6.7797, df = 3, p-value = 0.07926
##
## Model df: 6. Total lags used: 9The residuals are like white noise.
e. Use your chosen ARIMA model to forecast the next four years.
fc_gdp.NA_arima <- forecast(
gdp.NA_arima, h = 4
)
autoplot(fc_gdp.NA_arima)
The forecast values of gdp production are first increasing. But then
increase will stop and decrease a bit, before dampening to constant.
f. Now try to identify an appropriate ETS model.
gdp.NA_ets <- ets(gdp.NA)
gdp.NA_ets## ETS(A,N,N)
##
## Call:
## ets(y = gdp.NA)
##
## Smoothing parameters:
## alpha = 0.0586
##
## Initial states:
## l = 0.017
##
## sigma: 0.0129
##
## AIC AICc BIC
## -890.9744 -890.8931 -879.8731Chosen model is ETS(A, N, N).
g. Do residual diagnostic checking of your ETS model. Are the residuals white noise?
checkresiduals(gdp.NA_ets)
##
## Ljung-Box test
##
## data: Residuals from ETS(A,N,N)
## Q* = 32.364, df = 8, p-value = 8.014e-05
##
## Model df: 0. Total lags used: 8The residuals are not like white noise.
h. Use your chosen ETS model to forecast the next four years.
fc_gdp.NA_ets <- forecast(
gdp.NA_ets, h = 4
)
autoplot(fc_gdp.NA_ets)
The forecast values of difference in growth rate in GDP is
stationary.
i. Which of the two models do you prefer?
I prefer ARIMA model. Because it has better estimate given the residuals behaving like white noise. Also, I think at the current state, with the US government trying to reduce inflation, there should be slight decrease in the growth rate of GDP.
Question 9
beer=read.csv("beer.csv",header=T,dec=",",sep=";")
beer=ts(beer[,1],start=1956,freq=12) #,1 to take the values
lbeer <-log(beer) #make it has constant variance
train <- window(lbeer, end=c(1988,1))
h <- length(lbeer) - length(train)
ETS <- forecast(ets(train), h=h)
ARIMA <- forecast(auto.arima(train, lambda=0, biasadj=TRUE),
h=h)
STL <- stlf(train, lambda=0, h=h, biasadj=TRUE)
NNAR <- forecast(nnetar(train), h=h)
TBATS <- forecast(tbats(train, biasadj=TRUE), h=h)
Combination <- (ETS[["mean"]] + ARIMA[["mean"]] +
STL[["mean"]] + NNAR[["mean"]] + TBATS[["mean"]])/5
autoplot(lbeer) +
autolayer(ETS, series="ETS", PI=FALSE) +
autolayer(ARIMA, series="ARIMA", PI=FALSE) +
autolayer(STL, series="STL", PI=FALSE) +
autolayer(NNAR, series="NNAR", PI=FALSE) +
autolayer(TBATS, series="TBATS", PI=FALSE) +
autolayer(Combination, series="Combination") +
xlab("Year") + ylab("$ billion") +
ggtitle("log of Montly production of beer in Autralia")
c(ETS = accuracy(ETS, lbeer)["Test set","RMSE"],
ARIMA = accuracy(ARIMA, lbeer)["Test set","RMSE"],
`STL-ETS` = accuracy(STL, lbeer)["Test set","RMSE"],
NNAR = accuracy(NNAR, lbeer)["Test set","RMSE"],
TBATS = accuracy(TBATS, lbeer)["Test set","RMSE"],
Combination = accuracy(Combination, lbeer)["Test set","RMSE"])## ETS ARIMA STL-ETS NNAR TBATS Combination
## 0.07881405 0.07829543 0.08169119 0.09785364 0.07817191 0.07668909All of these individual models have really high accuracy, but the combination approach is even better.